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FOREWORD 


This  document  provides  a  detailed  description  of  the  formulation  and  computational 
processes  contained  within  the  Kozai  Mean  Element  Converter  (KMEC)  software  package. 
The  KMEC  program  was  developed  by  the  Naval  Surface  Weapons  Center  under  the  auspices 
of  the  Defense  Mapping  Agency  to  specifically  provide  operational  support  to  the  new  MX 
1502-DS  and  modified  TRANET  II  Doppler  beacon  satellite  tracking  receivers.  This  report 
has  been  reviewed  and  approved  by  Dr.  R.  J.  Anderle. 

Released  by: 

T.  A.  CLARE,  Head 
Strategic  Systems  Department 
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INTRODUCTION 


The  primary  function  of  the  Kozai  mean  element  converter  (KMEC)  is  to  generate 
mean  Kozai  orbital  element  sets  for  satellites  tracked  by  the  MX  1502-DS  and  modified 
TRANET  II  geoceivers.  These  geoceivers  will  use  the  mean  element  sets  in  conjunction 
with  a  Kalman  filter  to  provide  autonomous  orbit  updates  and  to  predict  station-satellite 
inview  periods.  To  avoid  problems  associated  with  mathematical  singularities,  such  as  those 
that  occur  with  near-circular  orbits,  the  following  nonsingular  element  set  has  been  selected 
to  perform  the  mean  element  transformation  in  KMEC: 


a  =  semimajor  axis 
A  =  S.  +  g  +  h 

£  =  e  cos  03  (co  =  g  +  h) 


where  a,  e,  i,  fi,  g,  and  h  are  the  usual  Keplerian  elements. 


(1) 


KMEC  is  comprised  of  seven  basic  computational  functions: 


1.  The  process  flow  supervisor  (PFS) 

2.  The  Cartesian  input  section  (CIS) 

3.  The  Brouwer  input  section  (BIS) 

4.  The  Walter  mean  element  interator  (WMI) 

5.  The  nonsingular  orbital  element  builder  (NEB) 

6.  The  geoceiver  format  section  (GFS) 

7.  The  Keplerian  element  builder  (KEB) 

Figure  1  shows  the  functional  overview  of  KMEC  and  the  interfunctional  data  flow.  A 
detailed  description  of  each  of  these  functions  is  presented  in  the  following  sections. 
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FIGURE  1.  KMEC  FUNCTIONAL  OVERVIEW  AND  DATA  FLOW 
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As  user-supplied  input,  KMEC  requires  the  selection  of  a  processing  option  (IOAOP). 
If  IOAOP  =  0,  KMEC  will  use  a  nonsingular  transformation  theory  to  convert  a  Brouwer 
mean  element  set  obtained  from  Space  Surveillance  (SPASUR)  input  data  to  a  Kozai  mean 
element  set.  If  IOAOP  =  1,  KMEC  will  use  the  same  nonsingular  transformation  theory  to 
convert  osculating  Cartesian  position  and  velocity  vectors  to  an  associated  Kozai  mean 
element  set.  Also  required  as  input  are  the  satellite  number  (KSA),  the  satellite  type 
(ITYPE),  the  satellite  frequency  offset  (AO,  and  the  universal  time  of  transit  of  the  vernal 
equinox  (tVE).  These  data  are  included  on  the  identifier  data  card.  The  Kozai  mean 
elements  and  associated  information  are  written  to  hard  copy  during  the  computational 
cycle. 


THE  PROCESS  FLOW  SUPERVISOR  (PFS) 

FUNCTIONAL  DESCRIPTION 

The  principal  tests  performed  by  the  PFS  are  to  receive  input  data,  direct  the  pro¬ 
cessing  flow,  and  output  computed  results.  Specifically,  the  PFS: 

1.  Receives  from  input  the  user-selected  processing  option 

2.  Receives  from  input  the  identifier  card  data 

3.  Receives  from  input  Brouwer  mean  elements  from  five-card  SPASUR  data  when 
IOAOP  =  0 

4.  Directs  processing  through  the  CIS,  BIS,  WMI,  NEB,  KEB,  and  GFS  functions 

5.  Converts  input  data  to  the  proper  computational  units 

6.  Writes  to  hard  copy  the  input,  intermediate,  and  output  mean  element  sets 

The  flow  of  the  PFS  function  is  presented  in  Figure  2. 

When  IOAOP  =  1,  the  CIS  function  is  entered  and  osculating  Cartesian  position  and 
velocity  vectors  are  received  along  with  the  vector  epoch.  These  vectors  are  converted  to  an 
osculating  Keplerian  element  set  and  are  used  to  initiate  the  Kozai  transformation  process. 


PROCESSING  EQUATIONS 

The  semimajor  axis  read  from  the  SPASUR  data  is  the  Kaula  semimajor  axis  aR 
expressed  in  earth  radii.  This  is  converted  in  the  PFS  to  the  Brouwer  mean  semimajor 
axis  a"  via  the  transformation 
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FIGURE  2.  PROCESS  FLOW  SUPERVISOR  LOGIC  FLOW 
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In  the  above  expressions,  a„  is  the  earth’s  semimajor  axis,  J„  is  a  zonal  harmonic  gravi- 

6  2 

tational  constant,  and  i"  and  e"  are  the  Brouwer  mean  inclination  and  eccentricity,  re¬ 
spectively. 


BROUWER  INPUT  SECTION  (BIS) 

FUNCTIONAL  DESCRIPTION 

The  Brouwer  input  section  accepts  the  Brouwer  mean  element  set  from  the  PFS 
and  converts  it  into  an  associated  osculating  element  set  at  epoch  tQ.  This  is  accomplished 
through  the  application  of  the  Brouwer-Lyddane  theory,1,2  which  has  been  modified  to 
include  the  effects  of  atmospheric  drag  (the  atmospheric  drag  decay  rates  are  nulled  during 
this  computation).  These  osculating  elements  are  then  used  to  formulate  the  osculating 
nonsingular  element  set  that  initializes  the  mean  element  interation  algorithm.  The  BIS 
processing  logic  flow  is  shown  in  Figure  3. 
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FIGURE  3.  BROUWER  INPUT  SECTION  PROCESS  FLOW 
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PROCESSING  EQUATIONS  FOR  THE  BROUWER-LYDDANE  METHOD 


The  equations  used  to  compute  osculating  orbital  elements  from  mean  Brouwer 
elements  and  associated  decay  rates  are  delineated  in  this  section.  First  define  the 


following: 


a"  = 

semimajor  axis  decay  rate 

e"  = 

eccentricity  decay  rate 

h  = 

time  rate  of  change  of  mean 

t 

time  from  epoch 

n  = 

O 

Gi/a"3)” 

V  = 

(1  -  e"’-)K 

0  = 

cos  i" 

72  = 

1/2  C20  a2/a"2 

I2'  = 

72JT4 

7z  =- 

C30  a2  a"~3r)~6 

74'  = 

-  3/8  C40  a4  a"“V8 

7s'  = 

-  C50  a2  a"- ^-10 

a  = 

1  -  502 

a  = 

1  -  1102  -  4O04  a'1 

7  = 

1  -  302  -  804  a"1 

6  = 

1  -  902  -  240 4  or1 

A  = 

1-5 02  -  1604  a-1 

(4) 


where  the  Cj0  (j  =  2,  3,  4,  5)  are  the  zonal  harmonic  gravitational  expansion  coefficients. 


Then  the  secular  terms  are  computed  from 


2"  =  n0t  |l  +  |  72'??(302  -  1)  +  —  y2'2V 


-  15  +  16r?  +  25i72 
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(5) 


+  (30  -  96r?  “  9Or?2)02  +  (105  +  144r?  +  25p2)04 
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g"  =  n(,t  2  y2>  “  +  32  y2'2  [~35  +  247?  +  257?2 


+  (90  -  192t?  -  126t?2)02  +  (385  +  360r?  +  45r?2)04] 


+  yj;  74*  [21  -  9t?2  +  (-270  -t  126t?2 )62  +  (385  -  189i?2)04]|  +  g/ 


and 


h"  =  n0t  373*0  +  y  72 '2  [(“  5  +  12t?  +  9r?2)0  +  (- 35  -  3 6r?  -  5t?2)03] 


+  -74'  (5  -  3t?2)0  (3  -  702)  >  +  h0" 


Tlie  long  period  (dependent  upon  g")  terms  are  computed  from 


=  35  7s.' 


8ie  =  ~  e"2  t?2  X  sin  i"  sin3g"  -  ~  (3y2'2  p  -  1074'  7)sin2  g" 

96  72  12  72 


_  P^T^r6"2  1,2  X  sin  sin  g"  +  a  ^  T3<  +  "i7  ys>  (4  +  3q'2^  5 

l/o  72  4  7  |_  16 


sin  i"  sin  g"  +  [372 '2  (3  -  IO7/  7 

2472  L  J 


9!  +  g'  =  g"  +  9."  +  >/2 


Hr!- 


372'2l2  +  e"2-  11(2  +  3e"2)02 


-  40(2  +  5e"2)04of 1  -  400e''2  06  or2 


+  1074  S2  +  e"2  -  3(2  +  3e"2)02  -  8(2  +  Se'^or1  -  8Oe"206oT 2 
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♦  *. 

72 


'2 


72  „  5  / 

— —  0  -  -  74  7 

4  6 


sin  2g"  +|  ~~  j]3e"  X  sin  i" 


35  7s' 


X  1-  e"(3  +  2e"2 )  sin  i"  + 


e"302 


1152  72'  L  I  '  sin  i" 

+  2e"302  sin  i"  i  5  +  3202  a"1  +  8O04  of2 


+  <- 


73e"02 


472 '  sin  i"  64 

\ 

(26  +  9e"2 ) 


Ill  L  e«  JL 

72 '  L  sini" 


cos  3g" 

(4  +  3e"2)  +  e"sini" 


8  ~  e"02  sini"  (4  +  3e"2) 

32  72 


(3  +  1602ar 1  +  4O04af2)  +  'A  ^7  sini". 

72 


3  -  e"2  (3  -  e"2) 


1  +  173 

e"(-  32  +  81e"4) 

[_4  +  3e"2  +  r?(4  +  9e"2)  J 


,  5  75'  2  o 

+  77  — 7  Vi 5 
64  72 


sin  i"  >  cos  g" 


(9) 


=  h"  +  — ■ ^-7— 7  —  <—  sin-1  i"  +  sini"  p  +  3202ar1  +  8O04a-2 
14472  1 2  !_ 
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11  +  800 2 a-1  +  2000 4 a” 2 
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+  IO74' 
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The  short  periodics  (dependent  upon  E ,  f  ,  £  )  are  computed  from 


a  =  k"t  +  a"  -  a"  (3 62  -  1)  + 

V 


3  v  7  (1  -  e"  cos  E')3 


[362  -  1  +  3  sin2 i"  cos(2g"  +  2f')] 


e  =  e"  +  e"t  +  8l&  +  ^  \-^—3  J3  -  e"2  (3  -  e"2) 

2  (  V6  1  +  V3  \ 

+  |s  +  e"  cos  f'  •  (3  +  e"  cos  f')j  cos  f'  +  — — - — — 


e"  +^3  +  e"  cos  f'  (3  +  e"  cos  f')j  cos  f '  cos  (2  f'  +  2g") 

-  r^LL-  (1  -  92)  [3  cos  (2g"  +  f')  +  cos(2g"  +  3f')l 
2  L  J 


*?272; 


i  =  i" - ; - r5ie  +  0  sin  i"  sin  f*  sin(2f'  +  2g") 

7j2sin  1 

3 

+  2&"y2Q  sin  i"  •  cos  f'  cos(2f'  +  2g")  +  -  h'd  sin  i"  cos(2f'  +  2g") 

g  +  £  =  g'  +  £'  +~~  j-  6a(f '  -  £"  +  e"  sin  f')  +  (3  -  502) 

^3  sin(2f'  +  2g")  +  3e"  sin(2g"  +  f')  +  e"  sin(2g"  +  3f')]j 

+  17TT  -  1)  (0  +  1)  sin  f'  +  3(1  -  02) 

4(1  +  v)  ( 

^(1  -  a)  sin  (2g"  +  f')  +  (a  +  1/3)  sin  (2g"  +  3f')jj 
h  =  h'  +  ^2e"72^  cos  f'  +  ^  72'0jsin(2g"  +  2f') 


-  e"y2'0  sin  f '  cos(2f'  +  2g")  -  372'0  (f'  -  £"  +  e"  sin  f') 
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and 


2  72 '  (4 


e5£  =  r  \  ~  72^  -  7  74  7f  sin  2g" 


-  |  7  ^  t?3  sin  i"  +  4 V3  sin  i"  (4  +  9e2")5  1 
(4  72'  64  72'  ) 


cos  g"  +  f?3e"  X  sin  if  cos  3g" 

384  72 ' 


(16) 


-  ~  72  V  |  2(302  -  1)  (a  +  1)  sin  f'  +  3(1  -  02)  (1  -  a)  sin  (2g"  +  f') 

+  (o  +  -)  sin  (2g"  +  3f ')  | 

3  i  ) 


where 


The  eccentric  anomaly  E'  is  obtained  from  a  Newton -Raphson  iteration  upon  the  Keple; 
equation 


E'  -  e"  sin  E'  =  fi" 


(is; 


and  the  true  anomaly  f'  is  found  from 


sin  f'  = 


r\  sin  E* 

1  -  q"  cos  E' 


(19 


cos  f' 


cos  E'  -  e" 

1  -  e"  cos  E' 


(20 


11 
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The  final  osculating  values  for  a,  i,  and  h  are  computed  from  Equations  11,  13,  and  15, 
respectively.  Equations  5,  12,  14,  and  16  are  used  to  calculate  final  osculating  values  for 
12,  g,  and  e  for  the  following  relations: 


A  =  e  cos  J2"  -  &8Z  sin  Z" 

(21) 

B  =  e  sin  Z"  +  e5£  cos  Z" 

(22) 

Z  =  tan”1  (B/A) 

(23) 

g  =  (fi  +  g)  -  Z 

(24) 

and 

e  =  (A2  +  B2)*  (25) 

CARTESIAN  INPUT  SECTION  (CIS) 

FUNCTIONAL  DESCRIPTION 

The  primary  tasks  performed  by  the  CIS  function  are  to 

1.  Receive  from  input  inertial  Cartesian  position  and  velocity  vectors  at  epoch  to, 
i.e.,  r(to)  and  r(to) 

2.  Transform  the  osculating  inertial  Cartesian  components  to  osculating  Keplerian 
orbital  elements 

The  process  flow  of  the  CIS  function  is  shown  in  Figure  4. 

PROCESSING  EQUATIONS 

The  osculating  inertial  Cartesian  vectors  r(tQ)  =  (x,  y,  z)  and  r(to)  =  (x,  y,  z)  are 
transformed  to  osculating  Keplerian  orbital  elements  by  using  the  following  relationships. 


THE  NONSINGULAR  ORBITAL  ELEMENT  BUILDER  (NEB) 


The  NEB  function  uses  the  osculating  Keplerian  elements  obtained  from  either  the 
BIS  or  CIS  functions  to  form  the  osculating  nonsingular  element  set  given  by  Equation  1. 
This  nonsingular  element  set  is  used  to  initialize  the  Walter  mean  element  iterator  discussed 
in  the  following  section.  The  NEB  process  flow  is  shown  in  Figure  5. 


{ .TW7"*.  *  * « *  *  ,  *  - ,  -  \  u  vrv*  v  „>  **.  .-v 

£s  i'«  -V*"  *v*  ; ",r  -  '  ‘ S’ -‘"'Vr 

V, V-  >>  ■  *-V-  V.%’‘ oV  V  4  ■,  <,-■  t.W  -/ 


V.vS 


•  »  »  ”< J v .  •  .  *  *  ■ 
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FIGURE  5.  THE  NEB  FUNCTION  PROCESS  FLOW 
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THE  WALTER  MEAN  ELEMENT  ITERATOR  FUNCTION  (WMI) 

FUNCTIONAL  DESCRIPTION 

The  primary  tasks  performed  by  the  WMI  are  to 

1.  Iteratively  solve  for  the  mean  nonsingular  element  set  associated  with  the  osculat¬ 
ing  nonsingular  element  set  formed  by  the  NEB  function 

2.  Write  to  hard  copy  a  message  describing  the  convergence  status  of  the  WMI 
algorithm 

The  mathematical  computations  performed  by  this  function  are  relatively  lengthy  and 
quite  complex.  They  are  described  in  detail  in  the  next  subsection.  The  WMI  process 
flow  is  depicted  in  Figure  6. 

PROCESSING  EQUATIONS 

The  iterative  technique  used  to  find  the  mean  nonsingular  element  set  from  the 
associated  osculating  elements  is  similar  t  •>  that  described  by  Walter.3  This  mean  nonsingular 
element  set,  represented  by  j3.,  where 


is  obtained  from  the  iterative  process  executed  according  to  the  scheme 


ENTER  WMI  FUNCTION 


4,0=  /if 


EVALUATE: 


f  (  il?2opq  )  dX 

J  an ( 

EQS.  57  THRU'  62 


COMPUTE: 


EQS.  37  THRU'  42 


EVALUATE: 
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in 


Jk)  sp  /  Oc - 1 )  _< k -  1  )\ 

Pi  (to)  =  w  -  Pi  W.  •••>  W ).  a  =  1, 2,  ...,  6), 


until  the  condition 


_  (k)  (k-i) 

Pj(t0)  -  ?j(t0)  <  e,  (i  =  1,  2,  ....  6) 


is  satisfied  for  all  j.  in  the  above  expressions,  k  is  the  iteration  counter;  0°(to)  and  0.  (to) 
represent  the  osculating  values  and  the  short  periodic  variation  of  the  jth  element  at  epoch 


to;  and  E.  is  the  convergence  tolerance  for  the  jth  element. 


To  initiate  this  iterative  process,  it  is  assumed  that 


P.  (to)  =  /}?( to),  (j  =  1,  2,  ....  6). 


When  the  condition  in  Equation  34  is  met,  then 


0j(to)  =  Pj  (tc),  0  =  1,  2,  ...,  6) 


As  can  be  seen  from  Equation  33,  short  periodic  variations  for  the  nonsingular  element 
set  are  required.  These  can  be  obtained  by  integrating  the  associated  Lagrange  planetary 
equations  using  only  the  J2  zonal  harmonic  in  the  gravitational  disturbing  function: 


3  3 

1  -  -  sin2i  +  -  sin2i  cos  2  (a?  +  0 
2  2 


1  -  -  sin2ij  (1  -  e2)"3/2 


and 


« s /(?-)*•«  5/e-)* 


QS1> 


1 

2n2a27 


+ 


1 

2n2a27 


(42) 


In  the  above  expressions,  f  is  the  true  anomaly 


7  -  Vi  -e2 


(43) 


a(l  -  e2) 

1  +  e  cos  f 


(44) 


and  R20)x)  is  the  (pq),h  contribution  to  the  geopotential  disturbing  function  due  to  the 
J  zonal  harmonic  and  is  given  in  general  by:4 


R 


J(0n,n  Kill  \&omnn  (APm  cos  0Pmn„  +  BSm  Sin02mpq)  + 


2mpq  g+ j  "  Snip  2pq  J  2mpq  v  2m  2mpq  2m 

3 


(45) 


^empq(A2m  S'n®8mpq  “  C0S^2mpq^ 


wnere 


(  Cfim  ’  £'m  even 

I  -S2m,  £-m  odd 


(46) 


(  S2m>  even 

'  c2m,  2-m  odd 


(47) 


*  v  i  ■  -■  v-  Vi ' .*  ^r*  x\  jt j  c  >r;  -r*:  rc  *rv  <  * **•  JCTa  ^ 


R 


6 


v^  '• '  1  w"^  »"  *1  ^  jr^  *rr  jr'j  ?g* 

‘  . tv-1 
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and 


06mpq  =  (fi  -  2p  +  q)X  -  me 


(48) 


Here  9  is  the  Greenwich  sidereal  time.  Note  that  for  the  case  =  2,  m  =  0 


A2o  ~  C20  =  J 1 
1^2  0  =  Sjo  =  0 


(49) 

(50) 


and 


02Opq  =  (2  -  2p  +  q)X 


(51) 


The  J8mp  and  K8pq  functions  in  Equation  44  are  the  inclination  and  eccentricity 


functions  given  by 


-  <-»“ 


(e+m)! 


j2 


2KP!(«-P)!  ph 


£  <--iy 


j  Q22-a-  2  j  (  j  _  Q2  y+  (a  -  |a  |)/2 


(52) 


K(7)  =  H)lql  28  (l+T)~C-,ql 


«  lil+x  k  (-l)r  / 2P- 22 V- 2PU-2pV 


r+t 


Spq 


Y'  Y-'  y"'  -  1  > 
2-j  2-j  2-j  r!  t! 

k  =  0  r=  0  t=0 


(53) 


(l+7)t+t-k  d_7)k}  (for  q  >  0) 


and 


|ql+k 


(-1)* 


KW  =(-l)W  £  •£  E-7^1 

Kpq  k  =  0  r=o  7^0  r,U 


-2p  W2p-2g\/g-2p+qN 
Jql+k-ry  \  k-t  )\  2 


r+t 


(54) 


(l+7)r+,-k  (1-T)k,  (for  q  <  0) 

21 


,-r_  St  r,  •*. 


-v-v-v-v 


N;/4 


^-5 

-is 

,  «  J 

'V, 

,o0 


L 


;  *  ‘  * 


r«  • 


iv; 

’rVVy 

I'SA 


m 


r  «  '  ' 

M 


g.»/ 


-U  >■ 

K>'"' 


c>:i 

'  -v-1 


yv 


L>  •  • 


}’>  vj 

jp“' 


1  -Toj 

f  J*  *  .  i 


’••V- 

I  * 


L 


; 


ft® 

St 

Lw\ 


'■■'j® 


k>< ; 
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where 


C=  cos 


k  -  integral  part  of 


jj  =  max(0,  -  a) 


£-m 

2 


j2  =  min(2£  -  2p,  £  -  m) 


a  =  m  -  £  +  2p 


The  <Rfimpq  and  IIfimpq  functions  in  Equation  45  are  given  by 
k  /|q|\  /  N  \ 

«fimpq  =  T.  J]  H)n  +  U  5u  (  \  ^l-u  7?"  pl“  l-2n  +  u  Qan-u 

n-o  ",  \  u  /  \2n  -  u/  cso 


k'  U  2 


»«»p,  *  E  E  <-'> 


n  =  0  t 

u  =  u  j 


aW  \  ,  , 

j^lq|-u  pM-2n-1+u  Q2 n  +  1  - u 
2n+l-u  / 


where 


Iql  +  W 


,  k'  - 


iql  +  ia|  -  1 


=  max(0,  2n  -  |a|),  u2  =  min(2n,  |q|) 


=  max(0,  2n+l  -  |a|),  u'  =  min(2n+l,  |q|) 


=  1,  if  q,  a  are  both  positive  or  negative 


=  (-  l)u,  if  q  or  a  is  negative 
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The  integrals  appearing  in  Equations  38  through  42  are  given  by  the  following 


expressions: 


^2  0  p  R2pq  |  ^20p<l  COS02Opq  +  R20pq  S*n  02Opq  | 


f/3R2  Opq\  , 

J\3X  )  ja  \  a3  /  j20p  ^2p<>  1  <^2°P<1  COS02Opq  * 

f/9R20pq\  /  3  \  (<\  . 

J\3a  /  d  ^2-2p+q y  2  y  a 4y  J2°p  ^pa  l^20^  Sm  0 


20pq  R20pq  C0S  02Opq  | 


aR2oPq\  /K\  t  /_J_\  /9K2pq\ 

3?  /  dX  2(  a3  )  20  p  ( 2-2p+q)  \3£  )  ^20p<i  +  lql  ^Pfl^opq'. 


Sitl02Opq-  1  !I20pq  +  ^  K2pq  ^Opq  '  C0S  02Opq 


IH  dx  -  - j2  p?  v  fak)  Un  ^  - qK- n- 


n,ftnn  +  q  K2pq  <i20pq'  1  cose 


/  \  /  i  \  [(/9J20p\ 

"  "  M  a3  y2™  l  2-2p+qJ  |\3P  J  20pq 


+  l2P-2l  J20p  ^mpq  t  Sln  0 


+  |2p  -  2|  J  U2m'  cos  0 
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(  q  -  1  (q  >  0) 

q'  =  \  (64) 

l  q  +  1  (q  <  0) 


-1  (2p  -  2  >  0) 

+  1  (2p  -  2  <  0) 


(65) 


The  partial  derivatives  of  the  inclination  and  eccentricity  functions  are  given  by 


UC^K'-JJcr-m^C^  JCJ'a'  J  -T  "i r^T^-jC e>--' j?*  jT_  iT- 
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3C 


where 


?!*  .  (.„«  -gaaL  r  <-.>  28 "  2p)(  2P 

3C  1  '  2*p!(S-p)!  2-,  1  y  \  j  /\s  -  m  -  j 

i=j, 


S“-I“l  [(2fi-la|)S2j  -  (2j+a-M)S2*-2]  \ 


S  ”  sin  (  2 


|  Q2t-a-2j-  i 


3  £  7  07 


3K*r<.  t,  9KsPq 


where 


v2pq  _  (-fi-|q|) 

r  1+y  *pq 


°°  |q  |+  k  k 


+  (-1)1^1 2s (1+7)  ~iMql  V'  V  V  - 

Lj  £_J  L~!  r!t 


k  =  0  r=  0  t=  0 


2p-2fi  \  /  -2p\  /  £_2p+q\  r+ 1 


lql+k-r  A  k-t 


k(l+rXl-7)k“1],  (for  q  >  0) 


(1+7)' 


x+t-k- 1 


[  (r+t-  k)(  1 — 7)k  - 


•  •*■  •***•«  « 


■  » -  -  >  * .  r  * 

*.*,.%  *,  j  ■ 

v  * ±  rl~  *.» 


S's' 

„•  j  **  4  *  *  *  «  r_- 
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3K 


gpq  _  (-g-|q|) 
37  1+7 


|q  l+k 


K 


8pq 


+  (-l)lql  28  (l+7)-8-|ql  E  EE 


t±y 

r!t! 


k  =  0 


r=  0  t=0 


/  -2p  \  /2p-2£ 

\|q|+k-r/  \  k-t 


fi-2p+qV+t 
,  2  / 


(l+7)r+t-k-  1 


t(r+t-k)(l-7)k  - 


(73) 


k(l+y)(l-7)k- 1 1  >  (for  q  <  0) 


and 


dKgp(2P-g)  _  -2£+l 
37  7  K«p(2p-2) 


27- 2 8+2 


E 


'  £-1  \ 
2k+|2p-  £1  / 


(74) 


2—  2  k  —  1 2  p — fi  |  k(  1 — 72  )k_1 ,  (for  q  =  2p  -  £  and  p'  =  — - -  ) 


It  should  be  mentioned  that  even  if  convergence  is  not  achieved  (Equation  34),  KMEC 
assumes  that  the  final  mean  element  values  are  correct  and  continues  processing  with  them. 
This  is  done  since  near  convegence  may  have  occurred  and  the  resulting  mean  elements 
may  still  be  usable.  Messages  concerning  the  state  of  nonconvergence  are  generated  to 
alert  the  user. 


THE  KEPLERIAN  MEAN  ELEMENT  BUILDER  (KEB) 

FUNCTIONAL  DESCRIPTION 

The  KEB  function  decomposes  the  mean  nonsingular  element  set  obtained  from  the 
WMI  function  into  a  mean  Keplerian  element  set.  The  KEB  process  flow  is  shown  in 
Figure  7. 
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PROCESSING  EQUATIONS 

The  decomposition  of  the  mean  nonsingular  element  set  into  the  mean  Keplerian 
element  set  is  accomplished  through  application  of  the  following  relations: 


e  =  (I2  +  n2)Vl 

(75) 

h  =  tan"1  (Q/P) 

(76) 

g  =  tan"1  (riff)  -  h 

(77) 

£  =  X  -  (g  +  h) 

(78) 

f=  2  sin"1  [(P2  +  Q2)Vj] 

(79) 

Of  course,  no  decomposition  of  the  mean  semimajor  axis  a  is  needed. 


THE  GEOCEIVER  FORMAT  SECTION  (GFS) 


FUNCTIONAL  DESCRIPTION 

The  GFS  function  assembles  the  satellite  ID,  type,  and  mean  Keplerian  or.bital  elements; 
computes  an  earth-fixed  longitude  at  epoch  for  the  mean  right  ascension  of  the  ascending 
node;  and  converts  the  epoch  from  modified  Julian  days  to  year,  day,  and  minutes  of  day 
(GMT).  These  data  are  written  to  hard  copy.  The  GFS  process  flow  is  illustrated  in  Figure  8. 


where  to  is  the  epoch  in  modified  Julian  days,  y  is  the  year  expressed  as  an  integer,  and 
d  is  the  day  of  year.  The  right  ascension  of  the  Greenwich  meridian  L  at  epoch  t  is 
computed  from 
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L  =  or  (t  -  t  ) 

©  v  o  ve7 


(81) 


where  co0  is  the  rotation  rate  of  the  earth  and  tye  is  the  time  of  transit  of  the  vernal 
equinox  expressed  in  modified  Julian  days.  The  earth-fixed  longitude  A  of  the  ascending 
node  of  the  orbit  is  then  computed  by  using 


A  =  h  -  L 


(82) 


REFERENCES 


1.  D.  Brouwer,  “Solution  of  the  Problem  of  Artificial  Satellite  Theory  Without  Drag,” 
The  Astronomical  Journal,  Vol.  64,  No.  1274,  1959,  pp.  378—397. 

2.  R.  H.  Lyddane,  “Small  Eccentricities  or  Inclinations  in  the  Brouwer  Theory  of  Art¬ 
ificial  Satellites,”  The  Astronomical  Journal,  Vol.  68,  No.  8,  1963,  pp  555-558. 

3.  H.  G.  Walter,  “Conversion  of  Osculating  Elements  into  Mean  Elements,”  The  Astro¬ 
nomical  Journal,  Vol.  72,  No.  8,  1967,  pp.  994-997. 

4.  G.  E.  0.  Giacaglia,  “The  Equations  of  Motion  of  an  Artificial  Satellite  in  Nonsingular 
Variables,”  Celestial  Mechanics,  No.  15,  1977,  pp.  191—215. 


NSWC  TR  83-135 


APPENDIX 

COMPUTER  LISTING  OF  PROGRAM  KMEC 


Cl  o  o  o  o  non  o  o  o  ooo 


,  .  •- .  <2  \  I  s'/' 
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PI  «  3.141  59  2  65  3  59 
DEGRAO  =  PI  /  180, 

REAO  *,  IOAQP 

READ  *,  KSA,  ITYPE*  FT,  T¥£ 

ENTER  CIS  FUNCTION 

IF*  I0A0P.NE.  0  )  CALL  ORBAOJt  IOAOF,KSA,UJO,  IYL0, 10 AY , SEC) 

IF  (  IOAOP.NE.O!  GO  TO  40 
REAO  4,  KSAP,  IYLD,IOAV 

4  FORMAT ( 2X, 15, 40 X, II , 13) 

REAO  5  ,  UJ0,  ETU,  HO,  GO,  BES,  BI 

5  FORMAT  C8X,F14.8,5(1X,F8.4>  » 

REAO  10,  AO 

10  FORMAT  1//8X, Fll.S) 

BXI  *  BI’OEGRAD 
BH  =  GO’DEGRAO 
BO  c  HO’DEGRAO 
BAM  =  ETU* DEGRAO 

SEE  EQS.  (2)  -  (3) 

FN  =  1.  /  (  1.  -  ES*ES  I  **  1.5 
TA  »  SI N(  BXI  ) 

TB  =  TA  *  TA 

TE  *  1.  -  (  3.*TBI/2. 

TAO  *  ((  3.  *  R2) / <2.*A0*A0 )  )*TE*FN 

BA  *  (  AO* ((1.  ♦  2.*TAD)  /  (  1.  -  TAO  ))  **  0.6666666667)  *  AE 
SEC  =  0.00 
PRINT  15 
15  FORMAT  (1H1  ) 

ENTER  BIS  FUNCTION 

30  CALL  BRAUER  »BO,82,83,B4,«5,DT,BA,8ES,BXI,BAM,BH,30,CN2, A,ES,XI, 
1  AN,H, 0,A1,E1,RN1> 

PRINT  34 

34  FORMAT {/»56X»*8R0UHER*) 

PRINT  35,  KSA ,UJD, AO, BA ,BES ,81, 8X1, ETU, BAM  ,G8 ,BH,MQ ,B  0 

35  FORMAT  <39X,*HEAN  ORBITAL  ELEMENTS  FOR  SATELLITE  *,I5,*«* 

1//9X, ’EPOCH  (JULIAN  DAY  MINUS  2  ,400,000.  5)  ♦,  E22. 14/ 9X , 

2*SEMIMA JOR  AXIS’, 24X.E22.14,*  EARTH  RADII  *,E22.14, 

3*  KILOMETERS*/9X, ’ECCENTRICITY* ,26X,E22, 14/9X, 

4* INCLINATION*  ,27 X, E 22.14, ’  QEGREES*,6X,E22.14,*  RADIANS*/9X, 
5*MEAN  ANOMALY*, 26X.E22. 14,*  DEGREES’, 6X,E22. 14,*  RADIANS*/9X, 
6*ARGUMENT  OF  PERIGEE*, 19X,E 22.14,*  DEGREES*,6X,E22. 14, 

7*  RADIANS*/9X, ’RIGHT  ASCENSION  OF  THE  ASCEADING  NODE  *,E22.14, 
fl*  DEGREES*. 6X,E22. 14,*  RADIANS*/) 

BA  a  A 

ENTER  NEB  FUNCTION 
40  CALL  FORM 

ENTER  UNI  FUNCTION 


A-4 


ooooooo  no  o  o  o  o 


^  KW  *CT!m  •CT  V\  O  v:’  XV  -_-3  .’CtXVxVX^C 


Ai.Vi  *  <  «  m  VX  <  ».  *  _  «.  *,  *.W  -• 
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CALL  MEAN 

XXXL  *  XM<2)  /  OEGPA0 
PRINT  15 

PRINT  50.  KSA,UJD,XM<1> ,  X  XX  L  ,  XM  (3)  ,  XN  (4» ,XM<5> ,XM<6 ) 

50  FORMAT (26X. ’MEAN  NOHSINGUt AR  ORBITAL  ELEMENTS  FOR  SATELLITE  *,I5, 
1*  -  *,//31X, ’EPOCH  (JULIAN  OAV  MINUS  2,400,000.5  »  *, E22 . 14/31X, 
2*SEMIHA JOR  AX  IS*, 24X.E22. 14 ,♦  K IL0NETERSV31X, ’LAMBDA*, 32X,E 22. 14, 
21X, 

3*0EGREES*/3iX»*ZETA’,34X,£22.14/3iX,*ET A*,  35X,E22«  1V31X,’P’,37X, 
4  E22.14/31X,*Q*,37X,E22.14/> 

ENTER  KEB  FUNCTION 

CALL  DC MP3 S 
XA  *  BA  /  AE 
XXI  =  XI  /  DEGRAO 
XH  a  M  /  DEGRAO 
XO  *  0  /  DEGRAO 
XAM  =  AM  /  OEGRAO 
PRINT  60 

60  FORMAT  (///,57X,*K0ZAI’I 

PRINT  35,  KSA.U JD, X A, 3A,ES, XXI, XI,XAM,AM ,XH,N, X0,3 

ENTER  GFS  FUNCTION 

CALL  MX15I 2(KSA ,ITYPE ,FT, IYLO.IOAY ,6A,XXI,ES,XH, X0 , XA M, UJO.TVE, 

1  SEC, I0A0P) 

STOP 

ENO 

SUBROUTINE  FORM 


THIS  IS  THE  NEB  FUNCTION. 

OSCULATING  NONSINGULAR  ELEMENTS  ARE  F0RME0  FROM  THE  OSCULATING 
KEPLERIAN  ELEMENTS.  SEE  EOS.  («. 


COMMON  /  K0R8EL  /  8A,  ES,  XI,  M,  0,  AM 
COMKGN  /  NORBEL  /  A,  XL,  Z,  XN,  P,  C 
PI  =  3.14*59265359 
PI2  =  2.«PI 

A  s  BA 

XL  *  M  ♦  0  *  AN 
XL  *  AMOOf  XL.PI2  » 

US  a  M  •  O 

MB  *  AHQEH  HBVPI2  ) 

Z  *  ES  *  COS(  MB  > 

XN  =  ES  *  S IN (  MB  ) 

P  =  SINC0.  5*XI)  *  C0S(  O  i 
Q  *  SIN  (0. 5’XI)*SIN(  O  I 
PRINT  9 

9  FORMAT  1/.46X, ’OSCULATING  NONSINGULAR  ELEMENT  SET—*) 
PRINT  10,  A  ,  XL,  Z,  XN,  P,  Q 
10  FORMAT </,50X»’A  =*, G16.  1 C,* KM* ,/50X, 
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♦LAM8D A=*,G16.  10,*RAD*,/5flX, 
♦ZETA  »*,G16. 1C,/5QX, 

♦ETA  -♦. Gl&« 10»/50X, 

♦P  *♦, G16. 1 0,/5QX, 

♦Q  =»,G16.  10) 


RETURN 

END 

SUBROUTINE 


THIS  IS  THE  HMI  FUNCTION.  MEAN  NONSINGULAR  ELEMENTS  ARE  OBTAINED 
USING  THE  HALTER  ALGORITHM. 


DIMENSION  X0SC<6), XHS6),T0L ( 6* , 0ELX<6> , XSP <6) , XMN< 6  ) 
DIMENSION  SLM(b) 

COMMON  /  NORREL  /  A  ,  XL,  Z,  XN  ,  P,  Q 
COMMON  /  MNEL  /  Xi 
COMMON  /  I NTG  /  SUM,  G,  XMOT 
OATA  TOL  /  0. 01, 5*0. 00001  / 

XOSC(l)  a  A 
XOSCt 2)  =  XL 
XOSC  (3)  =  Z 
X0SC( 4)  =  XN 
K0SC<5)  =  P 
XOSC (6)  «  Q 
KOUNT  «  0 


INITIALIZE  PROCESS.  SEE  EQ.  (35) 

DO  10  I  *  1,6 

10  XM(I )  =  XOSCtI) 


EVALUATE  SUMS  OF  INTEGRALS  OF  GEOPOTENTIAL  DISTURBING  FUNCTION 
PARTIAL  DERIVATIVES.  SEE  EQS.  (57)- (74),  (37)  -  (42). 

20  CALL  EVI 


LOAO  SHORT  PERIODIC  ARRAY. 

KSP (  1  )  *  XSPA  (  XH  ) 

XSP(2)  =  X  SPL (  XH,  SUM,  G,  XMOT  ) 

XSPC3)  «  XSP Z*  XH,  SUM,  G,  XMOT  I 

XSP(4)  a  XSPXNt  XM,  SUM,  G,  XMOT  ) 
XSP(5 )  »  XSFP(  XM,  SUM,  G,  XMOT  > 

XSP<6)  =  XSPQ(  XH,  SUM,  G,  XMOT  ) 

KOUNT  *  KOUNT  ♦  1 


ITERATE  FOR  MEAN  ELEMENTS.  SEE  EQS.  (33)-(34). 

DO  50  J  «  1,  6 

XMN(J)  b  XOSC(J)  -  XSP(J) 

DELX(J)  =  ABS  (  XMN(J)  -  XH(J)) 

IF  (  OELXt  J)  ,L£.  TOL ( J ) >  GO  TO  50 
40  XM ( J)  *  XKN(J) 
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50  CONTINUE 

00  60  K  *  1.6 

IF (DELX (K) • GT.T0L( K ) .ANO. KOUNT. LT.100)  CO  TO  20 
60  CONTINUE 

70  IF  (KOUNT  .LT.  100  )  GO  TO  9 0 
WRITE  (6,100)  KOUNT 

100  FORMAT ( /.2  4X  ,*THE  KOZAI  MEAN  ELEMENT  CONVERSION  ALGORITHM  OIO  NOT 

1  CONVERGE  IN  *,I5,*  ITERATIONS.*) 

PRINT  101,  (T0L(I)«I*1«6),(0ELX(I),X>1,6) 

101  FORMAT  (/,5X»*TH£  OESIREO  TOLERANCES  WERE--*,/, €(2X ,E1 2.  4)/5X, 

2  *THE  FINAL  TOLERANCES  WERE~*,/,6I2X,E12.4) ) 

RETURN 

90  CONTINUE 
PRINT  91,  KOUNT 

91  FORMAT (/,30X»*THE  KOZAI  MEAN  ELEMENT  CONVERSION  ALGORITHM  CONVERGE 
10  IN  * ,15, *  ITERATIONS.*) 

RETURN 

ENG 

SUBROUTINE  EVI 


THIS  IS  PART  OF  THE  HMI  FUNCTION.  SUMMATICNS  OVEi  P  ANO  Q  OF 
THE  INTEGRALS  OF  THE  GEOPOTENTIAL  DISTURBING  FUNCTION  PARTIALS  ARE 
EVALUATED.  SEE  EQS.  (6.26)  -  (57) -(74),  (37)  -  (42). 


DIMENSION  XM(6),  SUM(6) 

COMMON  /  INTG  /  SUM,  G,  XMCT 

COMM  CN  /  CON  /  OEGRAO,  XMU,  XJ2,  AE 

COMMON  /  MNEL  /  XH 

COMM  CN  /  KORBEL  /  A A, ES ,XI , H,OM , AM 

XMOT  *  SORT (  XMU/  <  XM(1)**3I) 

L  »  2 
M  *  0 

00  10  I  e  1,  6 
SUH(I)  =  0.0 
10  CONTINUE 
CALL  DCMPOS 
AA  *  XMU) 

G  =  1.  -  ES*ES 
G  *  SQRT(G) 

A  =  >N(3) 

B  x  XH( 4) 

C  »  XM<  5) 

D  =  XM(  6) 

CM  »  -XJ2*(XHU*AE*AE/(AA**3)) 

CHI  a  -CM/AA 
DO  20  1  =  1,3 

IP  *  I  -  1 
IQT  =  2*IP  -  2 
IAL  *  M  -  L  ♦  2*IP 
DO  30  J  =  1,  5 
IQ  x  j  -  3 

IF(IQ.EQ.IQT)  GO  TO  30 
THT  =  «L-2*IP+IQ)*XM(2) 
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IQP  x  0 

IFIIQ.GT.O)  IQP  x  IQ  -  A 
IF  (  IQ.LT.  0  »  IQP  ■  IQ  ♦  1 
MP  =  0 

IF  (  IAl  .LT.  0)  HP  a  M  ♦  1 

IF  I  IAL  .  GT  .  0  >  HP  x  N  -  1 

CJ  «  XJLMP <L*H«IP* XII 
CK  x  XKLPQ(L,IP,IQ, G) 

DKOZ  o  -  (  XM(3)/GI*OKLPQCL,  IP,  10,  G) 

DKDN  *  -  (  XN<4»/GI*OKLPQIL,  IP,  IQ,  G> 

0 JOP  =  -2. *XN  (5)*0 JLMP( L, H, I P,X I) 

OJDQ  «  -2. *XH(6>*DJLHPCL.M,IP,XI> 

R  «  RLHPQIL»H,IP,IQ,A,B,C,D) 

RQP  x  RLMPQU.»N,lP,IQPvA,B,C,0> 

RHP  x  RLMPQ  (L«MP,IP,IQ,A,B,C*D) 

BI  x  BILHPQ (L,H ,IP, IQ,A ,B,C ,0) 

BIQP  *  BILHPQ  (L«H,IP,IQP, A, B,C, D) 

81 HP  x  8ILHPQ IL»MP» IP, IQ, A, B,C, 0) 

CT  x  C0S(T HTl 
ST  x  SIN(THT) 

SUH (1)  =  SUHIl)  ♦  CM*CJXCK* IR*CT*BI*ST> 

SUM (2)  «  SUMI2)  ♦  <3./<2.-2.«IPHC>)*Ct‘l*CJ*CK*IR*ST-BI*CT» 

SUMI3I  x  SUH(3)  ♦  ( l./(  2.  “c  «*IP+IQ> ) *CH*C  COKOZ*R#-I  ABSIIQ)  *CK 
1  *RQP»*ST-( 0K0Z*BI*IABS II Q)*GK*8IQF)*CT> 

SUM  (4  )  «  SUM  141  ♦  U./(2.-2.*IP  +  IQI>*CN*CJ*UOKON*R-IQ*CK*BIQP|*Sr 
1  -  IDK0N*BI  +IQ*CK*RQP)*CT» 

SUH15)  x  SUMI5H-<l./(2.-2.*miQ>»*CH*CKMtOJDP*RH  ABSI2*IP-2>* 

1  C J*RHP) *ST  -  IDJDP*8I+IABSI2*IP-2»*CJ*8IMP»*CT) 

SUM! 6)  x  SUHI6)  ♦  I l./l 2.-2 .*IP*I Q) )*CM*CK«tt DJDQ*R-I 2,*IP-2) * 

1  CJ*BIFP>*ST-10JDQ*BI+I2*IP-2)*CJ*RHP)*CT) 

30  CONTINUE 
20  CONTINUE 
RETUti  N 
ENO 

SUBROUTINE  DCHPOS 
C 
C 

C  THIS  IS  THE  KE8  FUNCTION.  MEAN  NONSINGULA f  ELEMENTS  ARE 
C  DECOMPOSED  INTO  HEAN  KEPLERIAN  ELEMENTS.  SEE  EQS.  (75) -(79). 

C 

c 

DIMENSION  XMI6I 

COHHCN  /  MNEL  /  XI 

COMHCN  /  KOfiBEL  /  A,ES,XI,H,OM,AM 

PI  =  3.14159265359 

PJ2  »  2.*PI 

ES  =  SQRTIXKt3>*XHI3»  ♦  XM(4>*XM(4») 

OH  x  ARTKQ(XH(6),XM(5)> 

HB  =  ARTNQ(XN/4).XM(3n 
H  *  MB  -  ON 

IF  (  K.LT.  0.0»  M  *  P12  ♦  H 

XI  x  2.*ASINISQRT(XMI5>»XNt5>«*M(6l«XN<6>>> 

A  x  XH(l) 

AH  =  XH  (2)  -  MB 

IF  AH.LT.  0.0)  AH  *  PI2  «•  AH 
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RETURN 

ENO 

SUBROUTINE  BRAUER(B0»B2»B3» B4,B5,DT,A2P,E2P,CI2P,GL62PfG0ZP,H02P, 
1CN2. A.CE.CI.CL.  G.H. A0T»  RNDjESO) 


THIS  IS  THE  BIS  FUNCTION.  A  DRAG  AUGMENTED  BROUWER  >  LYCOANE 
THEORY  IS  USED  TO  GENERATE  OSCULATING  KEPLERIAN  ELEMENTS  FROM 
MEAN  BROUWER  ELEMENTS.  SEE  EOS.  (4) -(25). 


1004  A2P2=A2P**2  8RLY0050 

A2P4*A2P2**2  BRLV0Q60 

CN0=S«RT(B0/CA2P2*A2P))  BRLY0070 

E2P2*E2P**2  BRLY0060 

ETA=«SCKTU.-E2P2>  BRLY0090 

SIKEI»S IN<  C I2PJ  BRLY0100 

THET  A«COS( CI2P)  BRLY0110 

THET  A2*THETA**2  BRLV0129 

TMETA4=THETA2*»2  BRLY0130 

THETA6«THETA4«THETA2  BRLY0140 

CJ2s-B2/(2.*B0*A2P2)  BRLY0150 

ETA2=ETA**2  BRLY0160 

£TA3=ETA2*ETA  BRLY0170 

ETA4=ETA2**2  BRLY0180 

CJ21P=CJ2/ETA4  BRLV0190 

CJ3lP«B3/<  B0*A2P2*A2P*ETA4*ETA2I  BRLV0200 

C J41P«  <  3 . ♦ B4  > /( 8. *B  0*A2P4*E YA4* ETA4»  BRLY0210 

CJ51P«B5/(B0*A2P4*A2P*ETA4*»2*ETA2I  BRLY0220 

FUN1=3.*THETA2-1.  BRLY0230 

FUN2*=1.-5.*THETA2  BRLY0240 

SINEI2=SINEI**2  BRLY0250 

Ai=A2P»CJ2*FUNl  BRLY0260 

A0«-At/ETA3  BRLY0270 

A2*3.*A2P*CJ2*SINEI2  BRLY0200 

F1JN5=1.-11.*THETA2-(48.*THETA4)/FUN2  BRLY0290 

FUN6=  -FUN1-(6.*THETA  4>/FUN2  BRLYD300 

FUN4=  THETA  2/5 INEI2  BRLY0310 

FUN22«FUN2**2  BRLY0320 

CJ21 F2=CJ2 1P**2  BRLYO330 

EOlPc ((E2P*ETA2)*(3.*CJ21P2*FUN5-10.*  BRLY034D 

1CJ41P«FUN6>  )/(24.*CJ21P»  BRLYS350 

E21P=-2.*E01P  BRLV0360 

E31P*((35.*CJ51P*E2P2*ETA2*SJNEJT*<FUN2-lie.*THETA4)/R)N2>>/(96.*  BRLY0370 
1CJ21P)  BRLY0360 

EilPs-.7  5*E3lPM<.25*ETA2«SINEI)*{CJ31P4.3125*CJSiP*i4.*3.*E2P2)*  BRLY0390 
1 <1.-9«*THETA2-(24,*THETA4)/FUN2)>  J/CJ21P  BRLY0400 

Cl 0*-<E2P* THETA) /IETA2*SIN£I >  BRLVQ41Q 

CI2«CJ21P*TH£TA*SINEI*1 .5  BRLY0420 

CU=E2P*CI  2*.  66666667  BRLY0430 

FUN7= (-.5*ETA3*CJ21P»/E2P  8RLY0440 

CL21P=<ETA3/CJ21P»M.25*CJ21P2*FUN5-.63333333*CJ41P*FUN6>  BRLY0450 

CL12P=CN0*(1.+1.5*C  J21P«ETA*FUNl«-.09375*CJ21P2*ETAt- t-15.  +  16.*ETA+  BRLY0460 
125.*ETA2+(30.-96.*ETA-90.*E  1A2)  *THETAZ*U05.U44.»ETA*25.»,ETA2)*  BRLY0470 
2THETA4)*.9375*CJ41P*ETA*E2P2*(3.-30.*THETA2+35.*THETA4>>  3RLY0480 
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CL22Ps.5«CNQ*CN2 

G21P= (l./(  24. *C J21P))*I-3,*CJ21P2*(2 •♦E2P2-11«*I2.  •■3»*E2P2)*THETA2BRLY0508 
i-40.*(2.+5.*E2P2)*THETA4/FUN2-400.*E2P2*THETA6/FUY22IHO.*CJ4lP»  8RLY051Q 
2(2.  BRLY0520 

3+E2P2-3. ♦( 2 •♦3.*E2P2)*THETA  2-0. *(2.45.»E2P2I*THETA4/FUN2-80.*E2P2*BRLV0530 
4THETA6/FUN22) )  BRLY0540 

G12P=CN0*{-1.5*CJ21P*FUN24-.09375*CJ21P2*(-35.424.*:ETA425.»ETA2+  8RLY0550 

l(90.-i92.»ETA-126.*ETA2>*rHETA2+(385.*360.*ETA*45,  *ET  A2> *THETA4) ♦  BRLYQ560 
2.3125»CJ41P*(2i.-9.*ETA2M-27«.  ♦i2e.*ETA2>*THETA2M38  5.-18«.*ETA2)BRLYQ579 


3*THETA4>)  8RLY058C 

H2=i.5*CJ21P»THETA  BRLY0590 

H3*-2.*H2  BRLV0600 

H1=.66666667»E2P*H2  BRLY0610 

H31Ps<<35.*CJ51P*E2P2*E2P*THETA)/(144.«CJ21P>)M.5/SINEI*CFUN2-C  BRLY062B 
U6.*THETA4)/FUN2>«-SINEI*C5.  ♦(32.*THETA2>/FUN2*80.*THETA4/FUN22) >  BRLV0630 
HllP»-.25*  BRLY0640 

1  H3 1P*(  I.  25*E2P*THETA)  /  (CJ21P*SIKEI))*  (CJ31P+  .3125*0 J51P*  BRLYQ650 

2(4.«.*E2P2J*(l.-9.*THETA2-(24.*THETA4)/FUH2)»1.8r5*CJ51P*SINE12*  BRLY0668 
3(4.*3.*E2P2>M3.m6.»THETA2>/FUN2m8.*THETA4)/FJN22))  '  BRLY0678 

H21P=(E2P2*THETA>/(i2.*CJ21P)*(-3.»CJ21P2*Cll.M80.*Tt€TA2)/FUN2*  8RLY0680 
1 (200 .*THET  A4I  /FUN22H-10  .*C  J41P»  (3.  *(16  .*THETA2)  /FUN 24-  (48  .*THETA4)/BRLY0698 


2FUN22)) 


BRLY0700 


H12P=CN0*THETA*(-3.*CJ21P*.375*CJ21P2*(-5.*12,*ETA*9.*ETA2*<-35.-  BRLY0718 


136. ♦ETA-5.  *ETA2I*THETA2)U.  25*C  J41P*(5.-3.*ETA21*(3  .-/.♦THETA2)  I 
AID*CJ51P/C J2 IP 
AI02=FUN2- (16.*THET  A4)/FUN2 
C1*35./364.*AI0*£TA3»E2P»SIN£I*AI02 
AID3=THETA  2/S1NEI 
AID4=THETA  24SINEI 
E2P3=E2P2*E2P 

0  2*35./ 1152.  ♦AIOM  (-E2P*SINEI»(3.*2.*E2P2H-E2P3*AI0  3)»AI024 
12.*EZP3*AI 04#(5  •  ♦(  32 •♦THETA  2)/F UN2-f(80 .♦THETA4) /FUN 22 )) 
C3=1.-9.*THETA2-(24.*THETA4>/FUN2 
AID5=C J31P/C J21P 

C4*.25»AID5»(-E2P*AID3>+5./64.*AIDM-E2P*AI03*(4.«-3.*E2P2>* 
1E?P*SINEI* ( 26.+  9.  *E  2P2) )»C3-15.  /32.*AI0*E2P*AI 04M  4.+  3.*E2P2>* 
2(3.  +  (16.nHETA2»/FUN2*(40.*THETA4>/FUN22» 

C5=E2P/f l.+ET  A3)* (3  .-E2P2*(  3.-E2P2)  I 

C6*CE2PM-32.t81.*(E2P2»E2P2)) ) /( (4.*3.*E2F2» »ETA* (4. ♦9c»E2P2)) 
C7=.25*AI05*SINEI*C5+5./64. *03* AI0*ETA2*SIhEI*C6 
C8*-.25*AI05»ETA3»SINEI-5./64.*AI0*ETA3*SINEI*(4.»9.*E2P2»*C3 
0510  T=OT 

CL2P=CL12P ♦0T4CL22P*0T**2<'CL82P  ♦  RNO»OT*OT 
CL2P= AHOO( CL2Pt 6.23318530717  96 > 

IF (CL2PI520. 530 »  530 
520  CL2P=CL2F+6. 2031853071796 
530  G2P«G12P*DT*GB.’P 
H2P=H12P*0T4H02P 
SIN£G*SIN(G2P) 

C8SZ NG=COS (G2P) 

01E=SINEG'(SINEG*(E31P*SINEG+E21P)4E11P»*EC1P 
H1P»((H31?*SINEG*H21P)*SINEG+H11P)*CQSING4H2P 
GPLP=G2P+CL2P+.  5*(CL21P+G21P)*S  IN  (  2.*G2P)  4(01*02)  *COS  (3.*  G2P) 
l*(C4-«-C7>*C0  SING 
CL1P=CL2P 


BRLY0720 

BRLY0730 

BRLY0740 

BRLY0750 

BRLY0768 

BRLV0770 

BRLV0780 

BRLY0790 

BRLY0800 

BRLY0810 

BRLY0820 

BRLY0838 

BRLY0850 

BRLY0860 

BRLY0878 

BRLY0880 

BRLY0890 

BRLV0918 


8RLY0940 

BRLY0950 

BRLY0960 

BRLY0970 

8RLY8988 

BRLV0990 


BRLY1020 
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U*CL2P 

100  OELTAU*  <U-E2P*SIN<U)-CL2P)/  ( l.-E2P*CQS  (U) ) 

U*U»D£LTAU 

IF< A8S( OELTAU) -1. E- 10)  2  00,  100.100 
200  U=U-{U-E2P*SIN(U)-CL2P)/<1.-£2P*C0S<U)> 

E«U 

SINE1P»SINIEI 
C0SE1P=C0S <E) 

G1P*  G2P 

AOlVR*l./< 1.-E2P*C0SE1P) 

SINF1P*A0IYR*ETA*SINE1P 
COSFiP*AOIVR* ICOSE1P-E2P) 

F1P*ARTNQ( SINF1P.CQSF1P) 

IF  <ABS (F1P-CL2P)- 3. 1415926535896) 220,210.210 
210  STOP 

220  FUN3*(1.+CN2*T)**. 66666667 
C0SFG*G0SC2.*  CG1P+F IP) ) 

SINFG«SIN<2.*IG1P*F1P)) 

ADIVR2*A0IVR**2 

A0IVR3=A0I0R**3 

CI*CI2P+CI0*01E*CI1*SINF1P*SINFGM2.*CI1*CCSF1P*CI2)*CQSFG 

FUN8=F1P-CL1P*E2P*SINF1P 

0018  H«H1P*(2.*H1*C0SF1P*H2)*SINFG-H1*SINF1P*C0JFG»H3*-UN8 
KFUN*H/6. 2631 853071796 
FUN9=KFUN 

H«H-FUN9*6. 2831853071796 
IF (H)8022. 8023,6023 

8022  H«H*6. 2831853071796 

8023  A- A2P/FUN3  + A0*(A1*A2*CQ SFGI  * AOIVR  ♦  AOT*DT 
A 1 06=  A0IVR2*ETA2*ADIYR 
AI07=SIN(2.»G2P+F1P) 

AI06=SIN(2.*G2P*3.*F1P) 

01=.2S*CJ21P*(6.*<5.*THETA2-1.  )*FUN8M3.-5.*THETA2)  *t  3.*SINFG* 
13v*E2P*AI07  *E2P*AID6  )) 

02*.25*CJ21P*  <2.*<3.*THETA2-1.)*<AI06*1,  )*SINF1P»3. * U.-THET  A2)* 
1  < (-A 106 ♦ 1. >*AID7MAI06*. 33333333)* AI06) ) 

AI D9=C0S(2 • *G2P+F1P) 

AI010*COS12.*G2P*3.*F1P) 

D3=-ETA2*.  5*CJ21P*<l.-THETA2)M3.*AI09f  AIOIO) 

ETA6I«l./( ETA3*ETA3) 

04=ET  A  El* ( C5  +COSF1P* (3 .+E2P*C0SF IF*  C3. +E2P*C0SF IP) )) 
O5=ETA6I*CE2P+C0SFlP**3.f-E2P*COSFlP*J3.»E2PJ,C0SFlP)  J) 

D6*  ETA2*CJ2*.5*C<3.*T«ETA2-1.)*04+3.*(1.-TH£TA2)*05*COSFG)*03 
GAL*GPLP*01*(E2P*ET A2)/ tl.+ETA) *02 

CE«(E2P«1.  )*(1.*CH2*0T) **.66e66666666667*l.*01E*06  ♦  ESD*OT 
£DL«.5*E2P*CL21P*SIN  <2.*G2P)*C8*COSIN6*E2P*Ct*CO>  I3.*G2P)- 
1  ETA3*Q2 
AIOl  4=SINI  CL2P) 

AID15»COSICL2P) 

ESL=CE  *AI  D 14 *E  OL*  AI  0 15 
ECL»CE  *AID15-£0L*AID14 


BRLY1030 

8RLY1040 

BRLY1050 

BRLY1070 

BRLY1060 

BRLY1090 

BRLY1100 

BRLY1110 

BRLY1120 

BRLY1130 

BRLY1140 

BRLY115Q 


BRLY1210 

BRLY1220 

BRLY1230 

BRLY1240 

BRLY1250 

BRLY1260 

BRLY1270 

BRLY1260 


BRLY1310 

BRLY1320 

BRLY1330 

BRLY1340 

BRLY1350 

BRLY1360 

BRLY1370 

BRLY1360 

BRLY1390 

BRLY141Q 

BRLV1420 

BRLY1430 

BRLY1440 


BRLY1470 

BRLY1480 

BRLY1490 

BRLY1500 

BRLY1510 


CE=SQRT (ECL*ECL*ESL*ESL) 
CL=ARTNO(ESL»ECL) 

G* GAL-CL 

G*AHOOCG, 6. 28 31853071796) 


BRLY1520 

BRLV1530 

BRLY1540 

BRLY1550 


o  o  o  o  o  o  o 
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IF(G)8024,8025,8025 

8024  G  *  G  ♦  6.  2831853071796 

8025  RETURN 
END 

SUBROUTINE  COROELIRV.HV.Q, TAUE, U> 


THIS  IS  PART  OF  THE  CIS  FUNCTION.  OSCULATING  CARTESIAN  POSITION 
ANO  VELOCITY  VECTORS  ARE  TRANSFORNEQ  INTO  OSCULATING  tCPLERI AN 
ELEMENTS  .  SEE  EQS  .  (26) -(31). 


OIHENSION  RV(6),HV(6) 

X«RV (1) 

Y=RV  (2) 

Z*RV (3) 

XD0T=RV<4) 

YD0T=RV<5) 

ZDOT=RV (6) 

R*SQRT<X*X»V*Y+Z*Z) 

VSQ»XOO’,>XDOT*YOOT*YOOT*ZOOT*ZQOT 

HX*Y*ZDOT-Z*YBOT 

HY«Z*XOOT*X*ZOOT 

HZ=X*YDOT-Y*XDOT 

HeSQRT(NX*HX»HY^HY»HZ*HZI 

AA«1. /  (2  ./R-VSQ/Q) 

ESINU=(X*XDOT*Y*YOOT*Z*ZOOT >/SQRT(Q*AA) 

ECOSUsR*VSQ/Q-l. 

ESQ3£SINU*ESINU*E20SU*EC0SU 

E«SQRT(ESQ) 

ROOT=SQRT ( 1 .-ESQ) 

ANGLEI*ARTNQCSQRT(HX*HX*HY*HY),HZ> 

IF (ANGLE I- T  AUE) 1,5,  5 

1  SGNH  2*HZ/A8S  (HZ  ) 

ANGLEI=1. 57 07963257 994* II. -S6NHZ> 

OHEGA*0. 

IF  (E  -TAUE) 2.3.3 

2  ANOMALcARTNQ(SGNHZ*Y.X> 

GO  TO  9 

3  PERIG*ARTNQ (SGNHZ*Y , X)-ARTN 0 (ROOT* ESINU.ECOSU- ESQ) 

IF (PERIG)4 .7, 7 

4  PERIG=P£RlG+6. 2831853071796 
GO  TO  7 

5  ONEGA  =  ARTNQ (HX.-HY) 

IF(E-TAUE) 8,6,6 

6  PERIG*ARTNQ(Z*H*(£COSU-ESQ)MX*HY~Y*HX)*ROOT*ESIND,Z*H*ROOT* 
1ESINU-1 X*M  Y-Y*HX> * (ECOSU-ES Q> ) 

7  U* ARTNQ (ES INU .ECOSU ) 

ANONAL=U-ESINU 

GO  TC  10 

8  ANONAL*AR.rNQCZ*H,Y*HX-X*HY> 

9  E  *  0» 

PERIG«0. 

10  CONTINUE 

C  TAUss- ANOHAL*SQRT  (AA*AA* AA/Q) /3600. 


03/04/83 


BRLY1560 

BRLY1570 


EL  0060 


-.i: "  x  •a^i^»aTueuauBu.TuxB»L»^vsEkiaaaa5nm^a  SA-O-'U**? 
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TAUe ANONAL 

HV (i)*AA 

HV ( 2)*E 

HV(3)*ANGLEI 

HV<4)*PERIG 

H¥<5>=OMEGA 

H¥(6)«TAU 

RETURN 

END 

SUBROUTINE  ORSAOJt  IOAOP,  KSA,  TOA  ,  IATR.  IDJDA,  AOJSE  ) 


THIS  IS  PART  OF  THE  CIS  FUNCTION. 


COMMON  /  KORBEL  /  A  01, E i» XI 1 ,G01, Hfil.FL 0 1 
DIMENSION  RV(  6)  ,  HV(6) 

IF  (  IOAOP  .  EQ  .  0  )  RETURN 
READ  *  ,  I  AYR,  AD JO  A,  AOJSE 
READ  *,  RV(1),  RV(2),  RV<3) 

READ  *  ,  R V (4)  ,  R'.  (5),  RV(6) 

TOA  =  17210  V,.  ♦  367*IAYR  -  <7*IAYR)/4  ♦  ADJOA  ♦  ( AQJ  SE /964  UQ. ) 

1  -  2*00001. B 


Q  s  398600.8 
TAUE  *  l.E-06 

CALL  COROEL  (  RV,  H V ,Q,T AU; , U  > 

AOl  =  HV( 1 ) 

El  *  HV (2) 

XII  *  HV (3 ) 

G01  =  HV(4 ) 

HOI  e  HV(5) 

FL01  *  HV( 6) 

PRINT  15 

15  FORMAT  <  1H1  ) 

25  PRINT  20,  IOAOP 

20  RHAT142X, ’PROCESSING  OPTION  *,I4,*  SELECTEO*,/) 

PRINT  10,  I  AYR,  AOJOA,  AOJSE 

10  FORMAT  (39X,*THE  POST  ORBIT  AOJUST  EPOCH  IS-  YEAR* ,15,*  OAV  *, 

1  F6.1,*  SEC  *,G16.10) 

IDJDA  «  -,DJOA 
PRINT  11 

11  FORMAT  C/,39X,*IN?UT  POST  ORBIT  AOJUST  CARTESIAN  VECTORS  ARE  — ♦) 
PRINT  12,  (  RV(I),  1*1,6  > 

12  FORMAT ( /»39X,*X  =*,G16.10,*  KM*/39X, 

1  *Y  b*,G16.10,*  KM*/39X» 

2  *Z  a*,G16.10 ,*  KM*/39X, 

3  *XOOTs*,G16.1Q,*  KM  PER  SECV39X, 

<*  *Y00Ts*,G16.1Q,*  KM  PER  SECV39X, 

5  *ZDOr=*,G16.10 K N  PER  SEC*) 

RETURN 


ENO 


SUBROUTINE  HX  1502( KSA , ITYPE ,FT, IY , 10,3 A, XXI, ES, XH,X 0, XAM , UJO «Ttf E« 
1  SEC,  10  /OP) 
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THIS  IS  THE  GFS  FUNCTION.  THE  FOLLOWING  SHOULO  BE  NOTED* 

(1)  THE  INPUT  VALUE  FOR  TVE  SHOULO  BE  PERIODICALLY  UPDATEO 
FROM  THE  NAUTICAL  ALMANAC. 

(2)  THE  YEAR  COHPUTATICN  IAYR  ASSUMES  THAT  ALL  VECTOR  EPOCHS 
ARE  BETWEEN  THE  YEARS  1980  AND  1989  .  EPOCH  YEARS  OTHER 
THAN  THESE  HILL  NECESSITATE  CHANGING  1988  IN  THE  IAYR 
EXPRESSION  TO  THE  APPROPRIATE  OECAOE. 


DATA  HE  /  7.292115656E-05  / 

PRINT  10 

10  FORMAT  (  1H1  ) 

PI  *  3.14159265359 
PI 2  *  2.*PI 
OEG  =  PI  /  188. 

IF  (  IOAOP  .EQ.  0  >  GO  TO  15 
IAYR  *  IY 
TH  *  SEC  /  60. 

GO  TO  18 

COMPUTE  MINUTES  OF  OAY  (GMT)  •  SEE  EQ.  (80). 

15  IAYR  *  1980  ♦  IY 

TS  ■  UJO  -  1721044.  -  367*IAYR  ♦  <7*IAYR)/4  -  10  *  2480001.0 
TM  *  (TS  *  86400.)  /  60* 

COMPUTE  EARTH  FIXEO  LONGITUDE  OF  THE  ASCENOING  NODE.  SEE  EQ.  (81) 

18  RAG  =  HE  *  (  UJO  -  TVE  )  *  86400. 

RAG  *  AHOD  <  RAG.  P12) 

RAG  *  RAG  /  OEG 
XOL  =  XO  -  RAG 

IF  (  XOL.LT  .  0. 80  )  XX  *  360.  ♦  XOL 
PRINT  20 

20  FORMA T(////»37X,*GE00ETIC  SATELLITE  ORBIT  PARAMETERS  FOR  THE  MX  15 
102-DS  GEOCEIVER*) 

PRINT  30.  KSA .ITYPE .IAYR. ID .TM. XAM.XH.ES.BA. XOL. XX I .FT 
30  FORMAT ( //, 57X »*SATELLITE  IDENTIFICATION*. 6X.I5/57X, 

1  *SATELLITE  TYPE*. 16X.I5/57X, 

2  ^ELEMENT  SET  EPOCH  (GMT)*,/7GX, 

3  *VEAR*»7X»I5/77X» 

4  *D  AY*  »7X» 1 5/76X, 

5  *  MIN*.5X,F7.2///44X, 

6  *NEAN  A  NOH ALY*  <  2  OX  .£22. 1  4  .*  0EGV44X, 

7  *ARGUHENT  OF  PERIGEE*. 13K.E22.14.*  0EG*/44X» 

8  *ECCENTR1CITY*«2QX.E22.  14/44X, 

9  *SEMI-H AJOR  AXIS*,17X,E22.14,*  KMV44X, 

A  *LONGIT UOE  OF  ASCENOING  N00E*,5X,E22.14»*  »ES*/44X, 

B  *INCL INATI0N*<21X.E22.14«*  DEG*/44X, 

C  *TRANSH ISSION  FREQUE NCY* , 10X.E22. 14,*  PPH>.') 

RETURN 

END 

FUNCTION  XSPA(XM) 


u  wj ~-u  rxv»v;?/;' 


[A  * « ' 
r  v  *  4  ** 
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COMPUTE  THE  SHORT  PERIODIC  VARIATION  OF  THE  SEMI-1AJOR  AXIS. 
SEE  EQ.  (3-7). 


DIMENSION  XHC6) 

COMM  CN  /  (CORBEL  /  A,  ES,  XI,  H,  ON,  AH 
COMMON  /  CON  /  OEGRAD,  XMU,  XJ2,  AE 

F  =  AM  ♦  2  .*ES*SIN  <  AM)*  C5./4.)  *£S*ES*SIN  I2.*AM)*C13»/  12.  I  *ES*ES* 
1  SINI3.*AM) 

P  «  XMU)  *  (  1.  -  ES*ES  ) 

RRc  P  /  <  1.  ♦  ES*C0SCF)) 

XSPA  «  XJ2*CAE*AE/XHC1> >*< I < XH( 1) /RR)**3 . )*(1 .-1.5* SINC XI)* 
i  SIN(XI)  *  1.  5*SIN<XI)*SIN<  XI)*C0S(2.*H*2.*F))-Ci.-1»  5*SINCXI>* 

1  SINIXI)  )*U.-ES*ES)**C-i. 5)) 

RETURN 

ENU 

FUNCTION  XSPL  <  XM,  SUM,  G,  XMOT  ) 


COMPUTE  THE  SHORT  PERIOOIC  VARIATION  OF  LAMBDA.  SEE  EQ.  (38). 


O/.HENSION  SUM  (6),  XH(6> 

XSPL  *  -  { 2.  /  <XMGT*XMOT*XM(  1)  )  )*SUHC2)* I G/C2.  *X HOT* XMOT* XM  (1 

1  1  *XM  ( 1)  )  )*  JXM(3)  *SUMI3)*XH<4)*SUMC4I)*C1./C2.*XM0T*XN0T*XM11>* 

2  XH(i) *G  )  )*  (XM(5)  *SUM(5)*XM(6)*SUH<6)  ) 

RETURN 

Ol  0 

FUNCTION  XSPZ  <  XM,  SUM,  G,  XMOT  ) 


COMPUTE  THE  SHORT  PERIOOIC  VARIATION  OF  XI.  SEE  EQ.  (39). 


DIMENSION  SUN  <6)  ,  XNC6) 

XSPZ  »  -<3/<XMOT*XH!OT*XMCi)  *XH  U)  ♦  U.  *C  )  )  )  *XM(  3)*Sl)MC  IS-  CG/CXHOT* 

1  X"(0.  *Xim>*XM(i)  ))*SUHC4)-fl./<2.*XH0T*XM0T*XNCl)*XMCl)*G))* 

2  XM(4>*<XH(5)*SUH<5)*XM16)*SUM<6>) 

RETURN 

END 

FUNCTION  KSPXN  (  XM,  SUM,  G,  XMOT  ) 


COMPUTE  )HE  SHORT  PERIODIC  VARIATION  OF  ETA.  SEE  EQ.  (40). 


m 


'M 

Iw* 

V-'.--'. 

o 

(W-  -«~ 


*5 


ftr« 


DIMENSION  SUM  16),  XH(6) 

XSPXffa-CG/  (XMOT*XMGT*XMC.t.>*XMiI  ) * (1.  *G ) )  )*XMC4) *SJH  <1 )♦ C G/C XMOT* 
X  XMOT* 

XHU)*XtHl))  )*SUMC3)  ♦  Cl ./(2,*XN0T*XNQT*XNl 1) *XN Cl )*G) >*XM<3 ) 
2  *(KM<5) *$UM«5)*XMC6)*SUH (6)) 

RETURN 

ENO 

FUNCTION  XSPP  (  XM,  SUM,  G*  XMOT  ) 


-  *  .  n 

c>V-] 


K  h  M  wj 

V/k  r, 

O  •*  '  hJ 


m 

kj) 


A-l  5 


NSWC  TR  83-135 


39.44.  26  03/04/83 


COMPUTE  THE  SHORT  PERIODIC  VARIATION  OF  F.  SEE  EQ.  (41). 


DIMENSION  SUM (6),  XH ( 6) 

Cl  *  l./(2 .*XHOT*XMOT*XH(l )*XH( 1)*G) 

XSPP  a  -C1*XM<5J*SUN<1>  -  9.5*C1*SUH<6)  ♦C1*XM  (£>)*(  XM  (4)  *SUM  <3>- 
1  XH ( 3 ) *SUH( 4) ) 

RETURN 

END 

FUNCTION  XSPQ  I  XM,  SUM*  C.  XMOT  ) 


COMPUTE  THE  SHORT  PERIODIC  VARIATION  OF  Q.  SEE  EQ.  (42). 


DIMENSION  SUM  161,  XM(6I 

Cl  a  1./(2.«XM0T*XM0T*XM<1>  *XM(1I*G) 

XSPQ  ■  -C1*XM  (6 1* SUM (11  ♦  0 . S*C1*SUM( 51  ♦  C1*XM<6>*(  XHC4)*SUM( 3)- 
1  XM(3)*SUH<4) ) 

RETURN 

END 

FUNCTION  XJLMP(L,N,P,I) 


EVALUATE  THE  INCLINATION  FUNCTION  3.  SEE  EQ.  (52). 


REAL  JLMP,  I 

INTEGER  F,  ALPHA,  AALPHA,  P2 
C  *  COS (0.5*1) 

S  a  SIN (0,5*1) 

ALPHA  e  M-L*2*P 
MALPHA  =  -ALPHA 
AALPHA  a  I  APS (ALPHA ) 

LI  *  L-M 
K  a  0  .5*L1 

DO  5  JJJ  =  1,10,2 

IF  ( LI* EQ, JJJ)  K  =  K  ♦  1 
5  CONTINUE 
L2  a  2*L-2*P 
P2  a  2*P 
Jl  =  0 
J2  a  LI 
L3  =  L»N 
L4  a  L-P 

FHULT  =  CF ACT (L3)/ I  FACT (P)*F ACT  (L4)*(2 .**L I))* ( (-1. )**K  > 

IF  (MALPHA. GT.J1)  Jl  a  MALPHA 
IF  (L2.LT.J2)  J2=t2 
IF  (J1.CT.J2)  PRINT  It*  Jl, J2 

10  FORM  AT  ( 1H0 , 20X»*PRQGR AH  TERMINATED  IN  FUNCTION  JLMP - Jl  GY  J2— 

A  Jl  AND  J2  =*,215) 

IF  (J1.GT.J2)  STOP 
JLMP  «  0.00 
Jl  a  Ji+1 


A-16 


Aj.'-V-  C^-VA.’-VA'-'V  ■“ 

-  -  -A*  *.*  1-  .1  ■  a.  K  ^  aLml.  tL. -  *  ^  _ /  .  . >  >  V  .»  '  J-  • 


W-'V 

k* 


E 


•;*S  i  V  ~  *  k^L%.i'%.'!  V'M  VT%'| 
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J2  =  J2*I 

00  20  JJ  »  J1,J2 
J  =  JJ-1 
L5  «  Li-J 

F  a  U-l.l  **J)*BIN0HCL2,  J ) *BIN0M (P2 ,L5) *FMULT 

JLMP  x  JLNP*F*(C**(2*L-ALPMA-2*J)  )*(S*4(ALPHA-AALPHA*2*J)  ) 

20  CONTINUE 

XJLMP  «  JLMP 

RETURN 

END 

FUNCTION  0JLNP(L,M,P.1> 

C. 

c 

C  EVALUATE  THE  DERIVATIVE  OF  THE  INCLINATION  FUNCTION  NIT H  RESPECT 
0  TO  C0S(.5*I)  .  SEE  EQ.  (68). 

C 

0 

REAL  I 

INTEGER  P,  ALPHA,  AALPHA,  P2 
C  «  COS <0.5*1 1 
S  =  SIN  (0,5  *1) 

ALPHA  =  H-L*2 *P 
HALPPA  «  -ALPHA 

« alpha  x  i abs (Alpha) 

LI  *  L“M 
X  =  0,5*Cl 

00  5  JJJ  =  1,11,2 

IF  (LI, EQ, JJJ)  K  =  K  ♦  1 
5  CONTINUE 
L2  =  2*L-2*P 
P2  =  2*  P 
Jl  s  0 
J2  x  LI 
L3  *  L+H 
L4  a  L-P 

FMULT  «  (FACT (L3)/(FACT (P) *FACT (L4)*(2,**L)))* (C-i, l**K) 

IF  (MALPHA.GT ,J1)  .11  ■  NALPHA 
IF  (L2.LT. J2)  J2=L2 
IF  (J1.GT.J2)  PRINT  10,  Jl,  J2 

10  FORHAT ( 1H0  ,20X, *PROGRAN  TERMINATED  IN  FUNCTION  OJLMP - Jl  GT  J2-- 

A  Jl  AND  J2  =*,215) 

IF  ( jisGT,  J2)  STOP 
DLMP  »  0.0 
Jl  a  JU1 
J2  *  J2*l 

OO  20  J J  *  J1.J2 
J  «  JJ-1 
L5  =  LI- J 

F  »  ({-l.)**J)*BIN0M(L2, J)*BINCM(P2,L5)*FMULT 

OLHP  =  OLNP  +  F*(C**(2*L-ALPHA-2*J-1))*(S** (ALPHA- AALPHA))* 

A  (  (2*L-  AALPHA)  ‘‘ (  S**{ 2*  J)  )  -  ( 2*J* ALPH A-AALPHA) * (S •  * (  2»J-2))  ) 

20  CONTINUE 

DJLHP  «  DLMP 

RETURN 

END 


>  *> 


I, 


\-r. 

$ 
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FUNCTION  XKLPQCLt Pt  QtGAHA) 


EVALUATE  THE  ECCENTRICITY  FUNCTION  K.  SEE  EQS.  (S3)  -  (54). 


INTEGER  P»Q,RtT,AQ,RR,  TT,PP,PU,PUP 
ITST  *  5?*P  -  L 
XLPQ  a  0.00 

IF  (Q.EQ.ITSTJ  GO  TO  100 

Lt  *  -2*P 

L2  «  2*P  -  2*L 

FI  *  0„5*<L  -  2*P  ♦  Q) 

AQ  3  IABS(Q) 

F2  »  1.  ♦  GAHA 
F3  *  1.  -  GAHA 

F4  3  < <-l.)**AQ)*{2.**L»*(F2»*(-L-AO»» 

DO  10  KK  «  It  3 
K  *  KK- 1 
KU  =  AQ  ♦  KK 

00  20  RR  =  1 . KU 

o  .  eg  .  j 

B0  30  TT  s  lf  KK 
T  «  TT-1 
L3  s  AQ  K  -  R 
L4  *  K  ••  T 

IFfQ.GE.OI  BIFAC=t(-l.)**R)*BIN0H(L2,L3)*BIN0HCLl,L4 
A  > 

IFIQ.LT.0)BIFAC3( (-1.}**T>*BIN0M(L1,L3>*9IN0M(L2,L4> 
TERM=F4*(BIFAC/(FACT<R) *FACT(T) ) ) • (FI** CR*T) ) * CF2** (R +T-K))  * (F3**K 
A  > 

XLPQ  3  XLPQ  ♦  TERN 
TST  =  A8S (TERN I 
TST1  *  0.61  *  XLPQ 
TST1  3  ABS  t  TST1 ) 

IF( ( T ST .LT .TST 1) • ANO. (TST . NE» 0. 0  0) )  GO  TO  40 
30  CONTINUE 

20  CONTINUE 

10  CONTINUE 

PRINT  11,  LtPtQt3AMA,TST,TSTl 

11  FORMAT  J1H  0,2X,*KLPG  010  NOT  CONVERGE-  L  P  0  GAHA  TST  TST1**53I5i 
A  3F15.7,//) 

GO  TO  48 

100  IPP  =  I  ABS  ( ITST) 

LG  3  L-l 
PP  *  (L  -  I PP 1/2 
PU  3  PP  -  1 
IF  (PU.LT. 0)  GO  TO  40 
PUP  =  PU  *  1 
DO  50  KK  =  ltPUP 
K  =  KK  -  i 
L7  *  2*K  ♦  IPP 

XLPQ  3  XLP0MGAHA**(l.-2.*Ln*BINCM«L6,L7>*BIN0H  <L7,  K»M2.*M-L 
A  7)IMU.-GAHA*GANA)'>*K) 

50  CONTINUE 
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40  XKLPQ  «  XLPQ 
RETURN 
END 

FUNCTION  0KLPQ<L,P,Q,GAHA) 


EVALUATE  THE  DERIVATIVE  OF  THE  ECCENTRICITY  FUNCTION  K  KITH  RESPECT 
TO  GANNA  <  ■  SORT  (1.  -  E*E>  >  .  SEE  EQS.  (72) -(74). 


INTEGER  P,Q,R,T,AQ,HR,TT,PP,PU,PUP 

ITST  a  2*P  -  L 

XL PQ  *  0.00 

XL  x  XKLP8  <L»P»Q»GAHAI 

IF  (Q.EQ.ITST)  GO  TO  100 

LI  *  -2*P 

L2  »  2*P  -  2*L 

FI  «  0. 5*( L  -  2*P  ♦  Q> 

AQ  a  I A8S ( Q ) 

F2  a  1.  ♦  GANA 
F3  a  1.  -  GANA 

F4  a  (<-i.  >**AQI*(2.**LJMF2**<-L-AQ»» 

XI  a  ((  -L-AQ  )/  F2)*X1 
DO  10  KK  =  1,  3 
K  a  KK- 1 
KU  a  AO  ♦  KK 

DO  20  RR  a  1 , KU 
R  =  RR  -  1 

CO  30  TT  *  1,KK 
T  =  TT-1 
L3  a  AQ  ♦  K  -  R 
L4  *  K  -  T 

IF(Q. GE. 0 )  BIFACe  < 1-1 . >**R)*8IN0N CL2» L3 »*BINON CL1,L4 
A  ) 

IF(Q.LT.O)BIFACaf t-l.»aaT»*BIN0N«Ll,L3)*aiN0N<L2,L4) 
TERH=FV*<B  I FAC/ (FACT <R>  *FAC  T  <T ) ) )* (Fl*» (RfT) )* (F2**  <R*T-K-1> I* 

A  <<F3«-*K>*(R+T-K»-K*F2*(F3**(K-1»>> 

XLPQ  *  XLPQ  ♦  TERN 
TST  =  ABS<TERN> 

TST1  a  «.Ol  *  XLPQ 
TST 1  =  ABS  I  TST1) 

IF(CTST.LT.TSTl).ANO.tTST.NE.O.«n»»  3Q  TO  40 
30  CONTINUE 

20  CONTINUE 

10  CONTINUE 

PRINT  11,  L,P,Q, GANA, TST, TST1 

11  FORNAT  UH0,3Xt*KuPQ  DIO  NOT  CONVERGE-  L  P  Q  GANA  TST  TSTla*,3I5, 

A  3F15. 7,//) 

GO  TO  40 

180  IPP  =  l ABS (IT ST} 

XI  =  ((  -2*U1UGANA)  *  XI 
L6  a  L-l 

PP  =  (L  -  IPP)/2 
PU  a  pp  -  1 
IF  <PU.LT. 0)  GO  TO  40 
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PUP  =  PU  ♦  1 
DO  50  KK  *  l.PUP 
K  «  KK  -  1 
L7  =  2*K  +  IPP 

XLPQ  *  XLPQ-2.MGANA**(2  .-2.*L>)*6INOM(L6,L7)*BIMOH(L7,K>* 
A  (2.**(-L7)>*K*((l.-GANA*GANA»*MK-il) 

50  CONTINUE 
40  OKLPQ  «  XLPQ  ♦  XI 
RETURN 
ENO 

FUNCTION  RLHPQ( L» N» P.Q. A. B»  C.O) 


EVALUATE  THE  R  FUNCTION.  SEE  EQ.  (SS). 


INTEGER  ALPHA .AALPHA, P.Q, AQ,U1,U2,L,UU1,UU2,UU 
ALPHA  «  M  -  L  ♦  2*P 
AALPHA  «IA BS 1  ALPHA) 

AQ  *  IABS(Q) 

R  *  0.00 

K  *  0.5*( AQ  *•  AALPHA) 

KK  =  K  ♦  1 
DEL  =  1.0 
IPRO  «  Q*ALPHA 
DO  10  NN  =  1,  KK 
N  »  NN  -  1 
11=  2*N  -  AALPHA 
L2«  2*N 
U1  «  0 

IF  (LI.  GT.  0)  U1  *  LI 
U2  J=  L2 

IF  (  AQ  .LT.  L2)  U2  =  AQ 
UU1  «  U1  *1 
UU2  —  U2  F  1 

IFIU1.GT .02)  PRINT  11.  Ul,  U2 

11  FORMAT  (1H  0.3X. ’LOWER  BOUND  G7  UPPER  IN  SUN  OVER  U  IN  FUNCTION  RLN 
APQ  -  Ul  U2  «*,2I6»//I 
IF  (U1.GT.U2I  STOP 

00  29  UU  *  UU1.  UU2 
U  *  UU  -  l 

IF  (IPRO. IT . 0)  OEL  «  <-1.0l**U 
N1  *  2*N  -  U 

R  «  R  ♦  ((-l.O>*MN*U)>*OEL*eiNQM(AQ,U)*BIN31  (AALPHA, Nl>* 

A  (A** (AQ-U))*  (B**U>*(C**(AALFHA-2*N*U)I*(0**(2*N-U>) 

20  CONTINUE 

10  CONTINUE 
RLMPC  »  R 
RETUSH 
ENO 

FUNCTION  8 ILHPQ (L.K.P.Q.A.B.C.O) 


EVALUATE  THE  I  FUNCTION.  SEE  EO.  (56)- 
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INTEGER  ALPHA.AALPHA. P. QtAQ.Ult  UP. U.UUltUU2»UU 
ALPHA  «  H  -  L  *  2*P 
AALPHA  *IABS (ALPHA! 

AQ  >  IABS( QJ 
R  =  0.00 

X  *  0.5*<AQ  ♦  AAL»hA  -  1) 

KK  «  K  ♦  1 
DEL  =1.0 
I PRO  «  Q*ALPHA 
00  10  NN  =  It  KK 
N  «  NN  -  1 

L  1  *  2*N  -  AALPHA  ♦  1 
L2*  2*N  ♦  1 
U1  «  0 

IF  <L1.GT.  0)  U1  «  LI 
U2  *  L2 

IF  (  AO  .LT.  L2)  U2  ■  AQ 
UU1  =  U1  +  1 
UU2  *  U2  ♦  1 

IF  <  U1  .GT.  U2)  R  =  0,  0 
IF  (  UL  .GT.  U2>  GO  TO  30 
00  20  UU  =  UU1,  UUP 
U  *  UU  -  1 

IF  tIPRO.LT.OI  OEL  «  <-1.0>4*U 
IF  {(  IPRO.EQ.Q) . AND. (Q.LT.O))  0£L*<-1. 0>**U 
IFUIPRC.EQ.O) .ANO.CALPHA.LT. 81)  OEL  ■  C-1.0>**U 
N1  =  2*N  -  U  ♦  1 

R  =  R*<<-1.0)*MNfUH))*QEL*BINCMiAQ.U>*6INOM«AALPHA,M1>* 

A  (A**<  AQ-UM  MB**U)  *<C**<  AALPHA-2*N»U-1>)  *  (D*M  2*N~U*1 »  > 

20  CONTINUE 

10  CONTINUE 
30  CONTINUE 
BILMPQ  «  R 
RETURN 
END 

FUNCTION  BINON(MtN) 


EVALUATE  A  BINOMIAL  EXPANSICN  COEFFICIENT  . 


IF  (N.LT.Q)  8XN0M=  0.00 
IF(N.LT.O)  RETURN 
IF(N.EQ.O)  BINOM  *  1.00 
IF  (N.EQ.8)  RETURN 
J  «  P 

IF  1M.LT.0  >  M*  N-M-l 
L  *  H  o  « 

IF  !  L  .  IT  .  0  )  BINOM  *  0.08 

IF  t  L  .  LT  .  0  >  RETURN 

BINOM  =  FACT(M)  /  <  FACT  INI  *FACT(L) i 

IF  <  J.LT.O)  BINOM  «  (<-l.«)»*K)  *  BINOM 

M  =  J 

RETURN 
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ENO 

FUNCTION  FACT (K) 


EVACUATE  A  FACTORIAL. 


IF  (K.LT. 0 )  FRINT  10, K 

10  FORMAT ( 1H0 , 20 X, *PROGR AH  HAS  TERMINATED  DUE  TO  FACTORIAL  OF  A  NEGAT 
AIVE  INTEGER— -Kb*,  I  5» 

IF  (K.LT.O)  STOP 
FACT  >  1.0 
IF  (K.EQ.O)  RETURN 
DO  20  I>  1  »K 
FACT  *  FACT  *  I 
20  CONT.NUE 
RETURN 
END 
6/7/«/9 
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